%% Table C2: Parameter Estimates, Education Intergenerational Transitions
% Code to compute the parameter estimates (alphas and  betas) of 
% the education intergenerational transitions

clear all
close all

filename = 'TransitionMatrix';
sheet = 'Avg_Schooling';
sheetU = '1988_UrbanArea';
sheetR = '1988_RuralArea';


filename = 'TransitionMatrixRFC';
sheet = 'Avg_Schooling';
sheetU = '1988_UrbanArea';
sheetC = '1988_CityProper';
sheetF = '1988_Slum'; 
sheetR = '1988_RuralArea';

E = xlsread(filename,sheet,'C4:E29'); %Avg years of schooling of each location (row: state, column: location(city,slum, rural))

e = zeros(26,3,5);
e(:,:,2) = xlsread(filename,sheet,'G4:I29'); %Midpoint of 1 to 4 years of schooling (the same for child and dad)
e(:,:,3) = xlsread(filename,sheet,'K4:M29'); %Midpoint of 5 to 8 years of schooling (the same for child and dad)
e(:,:,4) = xlsread(filename,sheet,'O4:Q29'); %Midpoint of 9 to 11 years of schooling (the same for child and dad)
e(:,:,5) = xlsread(filename,sheet,'S4:U29'); %Midpoint of 12 or + years of schooling (the same for child and dad)

Prob_C = zeros(5,5,26);
Prob_C(:,:,1) = xlsread(filename,sheetC,'K96:O100'); %MA: Transition matrix of urban area
Prob_C(:,:,2) = xlsread(filename,sheetC,'K109:O113'); %PI: Transition matrix of urban area
Prob_C(:,:,3) = xlsread(filename,sheetC,'K122:O126'); %CE: Transition matrix of urban area
Prob_C(:,:,4) = xlsread(filename,sheetC,'K135:O139'); %RN: Transition matrix of urban area
Prob_C(:,:,5) = xlsread(filename,sheetC,'K148:O152'); %PB: Transition matrix of urban area
Prob_C(:,:,6) = xlsread(filename,sheetC,'K148:O152'); %PE: Transition matrix of urban area
Prob_C(:,:,7) = xlsread(filename,sheetC,'K174:O178'); %AL: Transition matrix of urban area
Prob_C(:,:,8) = xlsread(filename,sheetC,'K187:O191'); %SE: Transition matrix of urban area
Prob_C(:,:,9) = xlsread(filename,sheetC,'K200:O204'); %BA: Transition matrix of urban area
Prob_C(:,:,10) = xlsread(filename,sheetC,'K70:O74'); %MG: Transition matrix of urban area
Prob_C(:,:,11) = xlsread(filename,sheetC,'K83:O87'); %ES: Transition matrix of urban area
Prob_C(:,:,12) = xlsread(filename,sheetC,'K5:O9'); %RJ: Transition matrix of urban area
Prob_C(:,:,13) = xlsread(filename,sheetC,'K18:O22'); %SP: Transition matrix of urban area
Prob_C(:,:,14) = xlsread(filename,sheetC,'K31:O35'); %PR: Transition matrix of urban area
Prob_C(:,:,15) = xlsread(filename,sheetC,'K44:O48'); %SC: Transition matrix of urban area
Prob_C(:,:,16) = xlsread(filename,sheetC,'K57:O61'); %RS: Transition matrix of urban area
Prob_C(:,:,17) = xlsread(filename,sheetC,'K304:O308'); %MS: Transition matrix of urban area
Prob_C(:,:,18) = xlsread(filename,sheetC,'K317:O321'); %MT: Transition matrix of urban area
Prob_C(:,:,19) = xlsread(filename,sheetC,'K330:O334'); %GO: Transition matrix of urban area
Prob_C(:,:,20) = xlsread(filename,sheetC,'K213:O217'); %DF: Transition matrix of urban area
Prob_C(:,:,21) = xlsread(filename,sheetC,'K226:O230'); %RO: Transition matrix of urban area
Prob_C(:,:,22) = xlsread(filename,sheetC,'K239:O243'); %AC: Transition matrix of urban area
Prob_C(:,:,23) = xlsread(filename,sheetC,'K252:O256'); %AM: Transition matrix of urban area
Prob_C(:,:,24) = xlsread(filename,sheetC,'K265:O269'); %RR: Transition matrix of urban area
Prob_C(:,:,25) = xlsread(filename,sheetC,'K278:O282'); %PA: Transition matrix of urban area
Prob_C(:,:,26) = xlsread(filename,sheetC,'K291:O295'); %AP: Transition matrix of urban area

Prob_F = zeros(5,5,26);
Prob_F(:,:,1) = xlsread(filename,sheetF,'K96:O100'); %MA: Transition matrix of urban area
Prob_F(:,:,2) = xlsread(filename,sheetF,'K109:O113'); %PI: Transition matrix of urban area
Prob_F(:,:,3) = xlsread(filename,sheetF,'K122:O126'); %CE: Transition matrix of urban area
Prob_F(:,:,4) = xlsread(filename,sheetF,'K135:O139'); %RN: Transition matrix of urban area
Prob_F(:,:,5) = xlsread(filename,sheetF,'K148:O152'); %PB: Transition matrix of urban area
Prob_F(:,:,6) = xlsread(filename,sheetF,'K148:O152'); %PE: Transition matrix of urban area
Prob_F(:,:,7) = xlsread(filename,sheetF,'K174:O178'); %AL: Transition matrix of urban area
Prob_F(:,:,8) = xlsread(filename,sheetF,'K187:O191'); %SE: Transition matrix of urban area
Prob_F(:,:,9) = xlsread(filename,sheetF,'K200:O204'); %BA: Transition matrix of urban area
Prob_F(:,:,10) = xlsread(filename,sheetF,'K70:O74'); %MG: Transition matrix of urban area
Prob_F(:,:,11) = xlsread(filename,sheetF,'K83:O87'); %ES: Transition matrix of urban area
Prob_F(:,:,12) = xlsread(filename,sheetF,'K5:O9'); %RJ: Transition matrix of urban area
Prob_F(:,:,13) = xlsread(filename,sheetF,'K18:O22'); %SP: Transition matrix of urban area
Prob_F(:,:,14) = xlsread(filename,sheetF,'K31:O35'); %PR: Transition matrix of urban area
Prob_F(:,:,15) = xlsread(filename,sheetF,'K44:O48'); %SC: Transition matrix of urban area
Prob_F(:,:,16) = xlsread(filename,sheetF,'K57:O61'); %RS: Transition matrix of urban area
Prob_F(:,:,17) = xlsread(filename,sheetF,'K304:O308'); %MS: Transition matrix of urban area
Prob_F(:,:,18) = xlsread(filename,sheetF,'K317:O321'); %MT: Transition matrix of urban area
Prob_F(:,:,19) = xlsread(filename,sheetF,'K330:O334'); %GO: Transition matrix of urban area
Prob_F(:,:,20) = xlsread(filename,sheetF,'K213:O217'); %DF: Transition matrix of urban area
Prob_F(:,:,21) = xlsread(filename,sheetF,'K226:O230'); %RO: Transition matrix of urban area
Prob_F(:,:,22) = xlsread(filename,sheetF,'K239:O243'); %AC: Transition matrix of urban area
Prob_F(:,:,23) = xlsread(filename,sheetF,'K252:O256'); %AM: Transition matrix of urban area
Prob_F(:,:,24) = xlsread(filename,sheetF,'K265:O269'); %RR: Transition matrix of urban area
Prob_F(:,:,25) = xlsread(filename,sheetF,'K278:O282'); %PA: Transition matrix of urban area
Prob_F(:,:,26) = xlsread(filename,sheetF,'K291:O295'); %AP: Transition matrix of urban area


Prob_R = zeros(5,5,20);
Prob_R(:,:,1) = xlsread(filename,sheetR,'K96:O100'); %MA: Transition matrix of urban area
Prob_R(:,:,2) = xlsread(filename,sheetR,'K109:O113'); %PI: Transition matrix of urban area
Prob_R(:,:,3) = xlsread(filename,sheetR,'K122:O126'); %CE: Transition matrix of urban area
Prob_R(:,:,4) = xlsread(filename,sheetR,'K135:O139'); %RN: Transition matrix of urban area
Prob_R(:,:,5) = xlsread(filename,sheetR,'K148:O152'); %PB: Transition matrix of urban area
Prob_R(:,:,6) = xlsread(filename,sheetR,'K148:O152'); %PE: Transition matrix of urban area
Prob_R(:,:,7) = xlsread(filename,sheetR,'K174:O178'); %AL: Transition matrix of urban area
Prob_R(:,:,8) = xlsread(filename,sheetR,'K187:O191'); %SE: Transition matrix of urban area
Prob_R(:,:,9) = xlsread(filename,sheetR,'K200:O204'); %BA: Transition matrix of urban area
Prob_R(:,:,10) = xlsread(filename,sheetR,'K70:O74'); %MG: Transition matrix of urban area
Prob_R(:,:,11) = xlsread(filename,sheetR,'K83:O87'); %ES: Transition matrix of urban area
Prob_R(:,:,12) = xlsread(filename,sheetR,'K5:O9'); %RJ: Transition matrix of urban area
Prob_R(:,:,13) = xlsread(filename,sheetR,'K18:O22'); %SP: Transition matrix of urban area
Prob_R(:,:,14) = xlsread(filename,sheetR,'K31:O35'); %PR: Transition matrix of urban area
Prob_R(:,:,15) = xlsread(filename,sheetR,'K44:O48'); %SC: Transition matrix of urban area
Prob_R(:,:,16) = xlsread(filename,sheetR,'K57:O61'); %RS: Transition matrix of urban area
Prob_R(:,:,17) = xlsread(filename,sheetR,'K304:O308'); %MS: Transition matrix of urban area
Prob_R(:,:,18) = xlsread(filename,sheetR,'K317:O321'); %MT: Transition matrix of urban area
Prob_R(:,:,19) = xlsread(filename,sheetR,'K330:O334'); %GO: Transition matrix of urban area
Prob_R(:,:,20) = xlsread(filename,sheetR,'K213:O217'); %DF: Transition matrix of urban area


%% Cities, Slums and Rural Regions
for j=1:5

%% function f= EstimatingTransitionParametersFunction(params, ProbData, Edata, educ_parent,R_F_C: R=1, F=2, C=3)

fun0=@(params) EstimatingTransitionParametersFunctionRFC(params, Prob_C, E, j,3);
                options = optimset('TolFun',1e-10);
                options = optimset('TolX', 1e-10);
                options = optimset('MaxIter',50000);
                options = optimset('MaxFunEvals',50000);
               % options = optimset('Display','iter','PlotFcns',@optimplotfval); figure
               PARAMS0= [0,0,0,0,0,0,0,0];
               [paramsC(j,:), fval_search_C(j,:),exitflagC(j,:),output]   = fminsearch(fun0,PARAMS0,options);

fun1=@(params) EstimatingTransitionParametersFunctionRFC(params, Prob_F, E, j,2);
                options = optimset('TolFun',1e-10);
                options = optimset('TolX', 1e-10);
                options = optimset('MaxIter',50000);
                options = optimset('MaxFunEvals',50000);
               % options = optimset('Display','iter','PlotFcns',@optimplotfval); figure
               PARAMS0= [0,0,0,0,0,0,0,0];
               [paramsF(j,:), fval_search_F(j,:),exitflagF(j,:),output]   = fminsearch(fun1,PARAMS0,options);
               
fun2=@(params) EstimatingTransitionParametersFunctionRFC(params, Prob_R, E, j,1);
                options = optimset('TolFun',1e-10);
                options = optimset('TolX', 1e-10);
                options = optimset('MaxIter',50000);
                options = optimset('MaxFunEvals',50000);
               % options = optimset('Display','iter','PlotFcns',@optimplotfval); figure
                PARAMS0= [0,0,0,0,0,0,0,0];
               [paramsR(j,:), fval_search_R(j,:),exitflagR(j,:),output]   = fminsearch(fun2,PARAMS0,options);
end


%% Reporting the Results and Constructing the Resulting Probability Functions

%% City Regions
RorU     = 3;
ProbData = Prob_C; 
Edata    = E;

for j=1:5,
    a2 = paramsC(j,1);
    b2 = paramsC(j,2);
    a3 = paramsC(j,3);
    b3 = paramsC(j,4);
    a4 = paramsC(j,5);
    b4 = paramsC(j,6);
    a5 = paramsC(j,7);
    b5 = paramsC(j,8);

for k=1:26  
      %  The Model
                D(1,k) = 1;
                D(2,k) = exp( a2+b2*Edata(k,1) );
                D(3,k) = exp( a3+b3*Edata(k,1) );
                D(4,k) = exp( a4+b4*Edata(k,1) );
                D(5,k) = exp( a5+b5*Edata(k,1) );
       Prob_Model(:,k) = D(:,k)./sum(D(:,k));  

       % The Data:
          % First, need to select whethere the region has parents with e(j)
          if sum(ProbData(j,:,k))==1  
             Prob_Data(:,k)= ProbData(j,:,k)';
             IC(k)  = 1;
             fk(k)  = sum((Prob_Data(:,k)-Prob_Model(:,k)).^2);
             fTk(k) = sum((Prob_Data(:,k)).^2);
          else,
             IC(k)   = 0;
             fk(k)   = 0;
             fTk(k)  = 0;
          end
end

         ObsjC(j) = sum(IC(:));
         fjC(j)   = sum(fk)/ObsjC(j);
         fTjC(j)  = sum(fTk)/ObsjC(j);
         R2Cj(j)  = 1-fjC(j)/fTjC(j);

end

        R2C       = 1-sum(fjC)/sum(fTjC);
        
clear fk fTk a* b*


%% Slums Regions
RorU     = 2;
ProbData = Prob_F; 
Edata    = E;

for j=1:5,
    a2 = paramsF(j,1);
    b2 = paramsF(j,2);
    a3 = paramsF(j,3);
    b3 = paramsF(j,4);
    a4 = paramsF(j,5);
    b4 = paramsF(j,6);
    a5 = paramsF(j,7);
    b5 = paramsF(j,8);

for k=1:26  
      %  The Model
                D(1,k) = 1;
                D(2,k) = exp( a2+b2*Edata(k,2) );
                D(3,k) = exp( a3+b3*Edata(k,2) );
                D(4,k) = exp( a4+b4*Edata(k,2) );
                D(5,k) = exp( a5+b5*Edata(k,2) );
       Prob_Model(:,k) = D(:,k)./sum(D(:,k));  

       % The Data:
          % First, need to select whethere the region has parents with e(j)
          if sum(ProbData(j,:,k))==1  
             Prob_Data(:,k)= ProbData(j,:,k)';
             IF(k)  = 1;
             fk(k)  = sum((Prob_Data(:,k)-Prob_Model(:,k)).^2);
             fTk(k) = sum((Prob_Data(:,k)).^2);
          else,
             IF(k)   = 0;
             fk(k)   = 0;
             fTk(k)  = 0;
          end
end

         ObsjF(j) = sum(IF(:));
         fjF(j)   = sum(fk)/ObsjF(j);
         fTjF(j)  = sum(fTk)/ObsjF(j);
         R2Fj(j)  = 1-fjF(j)/fTjF(j);

end

        R2F       = 1-sum(fjF)/sum(fTjF);        

%% Rural Regions
RorU     = 1;
ProbData = Prob_R; 
Edata    = E;

for j=1:5,
        a2 = paramsR(j,1);
        b2 = paramsR(j,2);
        a3 = paramsR(j,3);
        b3 = paramsR(j,4);
        a4 = paramsR(j,5);
        b4 = paramsR(j,6);
        a5 = paramsR(j,7);
        b5 = paramsR(j,8);
for k=1:20  
      %  The Model
                D(1,k) = 1;
                D(2,k) = exp( a2+b2*Edata(k,3) );
                D(3,k) = exp( a3+b3*Edata(k,3) );
                D(4,k) = exp( a4+b4*Edata(k,3) );
                D(5,k) = exp( a5+b5*Edata(k,3) );
       Prob_Model(:,k) = D(:,k)./sum(D(:,k));  
       % The Data:
          % First, need to select whethere the region has parents with e(j)
          if sum(ProbData(j,:,k))==1  
             Prob_Data(:,k)= ProbData(j,:,k)';
             IR(k)  = 1;
             fk(k) = sum(( Prob_Data(:,k)-Prob_Model(:,k) ).^2);
             fTk(k) = sum((Prob_Data(:,k)).^2);
          else,
             IR(k)   = 0;
             fk(k)  = 0;
             fTk(k) = 0;
          end
end

 ObsjR(j) = sum(IR(:));
 fjR(j)   = sum(fk)/ObsjR(j);
 fTjR(j)  = sum(fTk)/ObsjR(j);

 R2Rj(j)=1-fjR(j)/fTjR(j);

end


R2R=1-sum(fjR)/sum(fTjR);

paramsR
R2Rj

paramsF
R2Fj

paramsC
R2Cj



%% Constructing the Implied Transition Functions
paramsC =  [    0.3752    0.0602   -0.2782    0.0603   -1.4588    0.1549    2.2017   -0.8731
                0.3452    0.2609   -0.5171    0.4215    0.6049    0.2326   -0.6450    0.3876
                0.5579    0.0578    2.6716   -0.2078    3.5317   -0.1805    1.8001    0.0899
               -0.0527    0.2080   -0.3673    0.3916   -0.7364    0.7135   -0.9283    0.7899
               -2.9930    1.3062   -1.5812    1.2301   -1.2127    1.3262   -1.7434    1.5541];            

paramsF =  [       -1.3926    0.3519   -4.3194    0.8462    1.7253   -2.2400    5.6594   -6.0559
                    0.3826    0.1879   -1.3421    0.4963   -0.8398    0.1813    0.7335   -0.9144
                    0.1229    0.1407   -0.2644    0.3590    0.1099    0.1168   -0.3958   -0.2198
                    0.3636    0.0828    0.8828    0.0198   -0.1174    0.4303   -0.3179    0.1135
                   -0.5405    1.0020   -0.7428    1.0815   -0.8126    1.2632   -0.5466    1.0847];
           
paramsU =  [   -1.3447    0.3221   -4.0053    0.6585   -6.5037    0.9764   -0.2165   -0.5487
                1.5671   -0.0174    0.7875    0.1079    1.9935   -0.1216    1.2530   -0.0760
                2.1359   -0.0939    1.4595    0.0846    3.2487   -0.0719    1.8383    0.1454
                1.0022    0.0308    2.5224   -0.0673    3.3225   -0.0551    2.8724    0.0567
                2.9370   -0.0408    4.1377   -0.0805    4.3052    0.0418    3.8181    0.2640];
            

paramsR =   [   -1.4031    0.5588   -3.7842    0.7317 -153.0621    0.4551 -343.9319    0.6787
                 0.2595    0.2825   -2.5330    1.0551    3.1706   -3.5829  -11.6147    2.9642
                -5.0708    1.8372   -2.3179    0.9077    0.1668   -0.1033   -1.9433    0.5911
                 2.4999   -0.0915    1.4746    0.1801    3.4392   -0.4036    1.8471    0.2782
                 0.0341   -0.0626    2.8190   -0.9117    3.5031   -1.0195   -0.9551    0.6864];




